On some critical aspects of numerical simulation of seismic wave propagation

Numerical simulation of seismic waves is essential for understanding field observations and developing inversion schemes. The wavefields generated by the most popular numerical methods such as finite differences and finite elements are often contaminated with grid dispersion and edge reflections. We will present a grid dispersion and stability analysis based on a generalized eigenvalue formulation. Our analysis reveals that, for a spectral FE of order 4 or greater, the dispersion is less than 0.2% at 4-5 nodes per wave length and the scheme is isotropic. The FD and classical FE require a larger sampling ratio to obtain results with the same level of accuracy. The staggered grid FD is an efficient scheme but the dispersion is angle dependent with larger values along the grid axis. On the other hand, spectral FE of order 4 or greater is isotropic with small dispersion making it attractive for simulations for long propagation times. Further, we will discuss a new time-space domain finite difference scheme with adaptive length spatial operators. We also employ a simple scheme to absorb reflections from the model boundaries in numerical solutions of wave equations. This scheme divides the computational domain into boundary, transition, and inner areas. The wavefields within the inner and boundary areas are computed by the wave equation and the one-way wave equation, respectively. The wavefields within the transition area are determined by a weighted combination of the wavefields computed by the wave equation and the one-way wave equation to obtain a smooth variation from the inner area to the boundary via the transition zone.